StemPanTox: A fast and wide-target drug assessment system for tailor-made safety evaluations using personalized iPS cells

Summary An alternative model that reliably predicts human-specific toxicity is necessary because the translatability of effects on animal models for human disease is limited to context. Previously, we developed a method that accurately predicts developmental toxicity based on the gene networks of undifferentiated human embryonic stem (ES) cells. Here, we advanced this method to predict adult toxicities of 24 chemicals in six categories (neurotoxins, cardiotoxins, hepatotoxins, two types of nephrotoxins, and non-genotoxic carcinogens) and achieved high predictability (AUC = 0.90–1.00) in all categories. Moreover, we screened for an induced pluripotent stem (iPS) cell line to predict the toxicities based on the gene networks of iPS cells using transfer learning of the gene networks of ES cells, and predicted toxicities in four categories (neurotoxins, hepatotoxins, glomerular nephrotoxins, and non-genotoxic carcinogens) with high performance (AUC = 0.82–0.99). This method holds promise for tailor-made safety evaluations using personalized iPS cells.


INTRODUCTION
To date, chemical toxicity studies have been primarily conducted by in vitro testing in cultured human cancer cell lines or in animals such as mouse, rat, and rabbit. However, because these tests are not on ''normal'' human cells, their applications are limited (Perel et al., 2007). In addition, the use of animals has become a major issue from the standpoint of animal welfare; in 2019, the U.S. Environmental Protection Agency (EPA) announced that research studies using mammals as well as funding for mammal studies will be cut by 30% by 2025 and abolished by 2035 (Grimm, 2019). To reduce animal studies, new approach methodologies (NAMs) (Topical Scientific Workshop -New Approach Methodologies in Regulatory Science, 2016), which are any technology, methodology, approach, or combination thereof that can be used to provide information on chemical hazard and risk assessment that avoids the use of intact animals (New Approach Methods Work Plan: Reducing Use of Animals in Chemical Testing, 2016), are widely investigated and adopted for the development of new tools for testing pharmaceuticals and other chemicals for potential adverse impacts on human health and ecological endpoints, under the direction of the U.S. EPA and the amended Toxic Substances Control Act (TSCA) for improving regulatory toxicology guidelines to consider the use of toxicity information, computational toxicology and bioinformatics, and high-throughput screening methods and prediction models, prior to animal studies (Becker, 2019).
The embryonic stem cell test (EST) reported by Scholz et al. was the first to examine embryotoxicity in vitro using mouse fibroblasts, embryonic stem (ES) cells, and cardiomyocytes differentiated from ES cells; these developmental toxicity tests were previously performed only in animals (Scholz et al., 1999;Seiler and Spielmann, 2011). Later, this method was approved as a scientifically valid alternative by the European Centre for the Validation of Alternative Methods (ECVAM) (https://tsar.jrc.ec.europa.eu/test-method/tm1999-01). However, the EST uses mouse cells, and species-specific differences must be clarified in order to extrapolate its results to human. Subsequently, another research group reported that in an EST based on a human cell system (hEST), homologous human and mouse neurodevelopmental gene expressions are similar (Schulpen et al., 2015). Moreover, there is an extensive collection of literature on the research and development of pluripotent stem cell models for predictive toxicology, including assays based on specific biomarkers such as Hand-1 and Sox17 (Kameoka et al., 2014;Suzuki et al., 2011), assays customized to specific toxicities such as teratogenicity (West et al., 2010), models for relevant pathways such as Wnt signaling or TGF-b (Kugler et al., 2015;Uibel et al., 2010), detected by various modalities such as luciferase reporter, metabolomics, high-throughput imaging (Kameoka et al., 2014;Kleinstreuer et al., 2011;Uibel et al., 2010), and others, which have been summarized in the review papers of Luz and Tokar (2018) and Kim et al. (Kim et al., 2019).
In 2012, the U.S. Defense Advanced Research Projects Agency (DARPA) and the NIH invested a huge amount of budget on a national project to promote the development of biomimetic systems, leading to rapid progress in the field (Low et al., 2021). These systems, which mimic the (adult) human body and are constructed by filling tissues created in individual compartments with culture fluid and connecting them together, are expected to be used in human toxicity testing systems as an alternative to animals. Currently, however, very little progress has been made in adapting these systems for developmental toxicity testing, and no evaluation method has been established to determine how accurately these systems mimic the function of the normal human body (Allwardt et al., 2020). Realizing the practical application of these systems as high-throughput toxicity screening tools is likely to take several years.
For many years, the prediction of chemical toxicity has been carried out using a method based on the physicochemical parameters of the chemicals, referred to as the quantitative structure-activity relationship (QSAR) (Schultz, 1999). Computational QSAR models have been widely used in various fields of predictive toxicology such as DILI severity in hepatotoxicity or hERG ligands in cardiotoxicity (Ekins, 2014). However, there is a limit to the predictive ability of QSAR. One reason is that the mechanism that induces toxic responses resides within the cell, so information about the chemical alone cannot predict the response. In this regard, it should be possible to detect toxicity more accurately by obtaining information about changes in gene expressions in the cells. In addition, as stem cells differentiate, only the genes essential for that lineage are expressed and conserved through DNA methylation (Reik et al., 2001), whereas in pluripotent stem cells, a very large number of genes are expressed, including transporters and transcription factors; consequently, pluripotent stem cells are superior to differentiated cells as a tool for comprehensively detecting toxic chemicals. In light of these considerations, we developed hEST-GN (human embryonic stem cell test with gene networks), a prediction method that uses information on feature gene networks based on massive gene expression datasets obtained by exposing human ES cells to toxic chemicals as input data for machine learning. Using this method, we achieved highly accurate predictions of developmental toxicity categories (Yamane et al., 2016).
In this study, we advanced the hEST-GN and achieved high performance predictions of 24 toxic chemicals in broad toxicity categories, including adult toxicity. Furthermore, by selecting induced pluripotent stem (iPS) cells, which can be used as an alternative to ES cells in toxicity testing, and using a gene expression database created from ES cells, we were able to develop a system called ''StemPanTox'' that can predict chemical toxicity using iPS cell data via transfer learning and ameliorate the ethical issues related to ES cells. With further improvement, StemPanTox will contribute to the development of tailor-made, individualized toxicity assessment/prevention using personalized iPS cells, which would have enormous clinical value.

Development of an hEST-GN library for 24 chemicals and prediction using iPS cells
A schematic of the chemical assay is shown in Figure 1A. The human ES cell line KhES-3 was exposed to a total of 24 chemicals in six toxicity categories [neurotoxins, hepatotoxins, cardiotoxins, two types of nephrotoxins (glomerular nephrotoxins and tubular nephrotoxins), and non-genotoxic carcinogens] at six concentrations, including vehicle (solvent) alone. The chemicals were carefully assigned to the toxicity categories by referring to previous reports (Table 1). Gene expression data were obtained by RNA-seq at two time points, 24 and 48 h, after exposure. At each time point, a principal component analysis (PCA) was performed using transcription factor genes, and a total of 20 genes from the top five PCs were extracted as feature genes. Using these genes, gene network libraries were created for each of the 24 chemicals using the graphical Gaussian model (GGM) (Goodal, 1991). Similarly, the screened iPS cells were subjected to RT-qPCR to obtain gene expression data for the same 20 genes and create gene network libraries.
Next, using ES cell library labels, a chemical toxicity prediction system trained by both libraries from ES and  (Vapnik, 2006) using support vector machines (SVMs) (Cortes and Vapnik, 1995).

Gene expression response database of ES cells for 24 chemicals
To obtain the largest amount of data about the expression of genes that were perturbed by the 24 chemicals, ES cells need to be exposed to chemicals at the maximum concentration that does not cause an excessive degree of cell death. To this end, we first performed ATP assays and then plotted regression curves to calculate inhibitory concentrations (ICs) by carrying out serial dilutions of stock solutions containing the maximum soluble concentrations of the 24 chemicals. ICs in the range of 0.1%-50%, at which cell death begins to be observed, were set as the maximum exposure concentration (Table S1, Figure S1). The estimated IC50 and 95% confidence interval (CI) for each chemical are shown in Figure 1B (Table S2). Serial dilutions were carried out, with the maximum exposure concentration set as 1/1 to obtain ½, 1/4, 1 / 8 , and 1/16 dilutions, and a six-step exposure including vehicle alone was performed and repeated twice, yielding a total of 6 3 2 = 12 samples for each chemical. We collected RNA 24 and 48 h after exposure to the 24 chemicals, performed transcriptome analysis, and generated gene expression datasets for a total of 12 3 24 3 2 = 576 samples.
To examine the characteristics of the 24 chemicals at the level of differentially expressed genes (DEGs), we selected transcription factor-related genes (GO: 0,006,351) from the 576 datasets (4,032 genes). After lognormalization, batch effect elimination, and repeat merging, we generated DEG sets for which differences between each exposure data and vehicle values were significant (FDR <0.01 and log 2 |FC| > 1) and presented them in a heatmap ( Figures 1C and S2). According to this analysis, the number of DEGs was higher at 48 h than at 24 h for all concentrations, and over time, more genes were up-or downregulated due to exposure to the chemicals. At both 24 and 48 h, valproic acid, a strong neurodevelopmental toxicant, elicited gene expression patterns that were clearly distinct from those of the other chemicals. Similarly, lithocholic acid, a mammalian bile acid and well-known carcinogen, yielded distinct expression patterns.

Construction of gene networks of the 24 chemicals by the graphical Gaussian model (GGM)
To obtain feature genes used in the prediction, we performed PCA on the basis of the exposure data of each transcription factor gene, expressed as a log fold change (LFC) in expression relative to vehicle after log-normalization and batch effect elimination. There were 3,200 genes for the 24-h samples and 3,255 genes for the 48-h samples. For both exposure times, it was difficult to clearly separate the chemicals by the toxicity categories using two-dimensional PCA ( Figures 1D and S3). Accordingly, we selected two genes with maximum positive and negative loading values, which were considered to contribute the most to the first to fifth PCs; at each time point, 20 genes were selected as feature genes (Table S3). Among the 20 selected genes in each group, only ACTR3 (Welch et al., 1997), which has been implicated in cell shape and motility, was common.
Using the 20 selected genes, we estimated sparse gene networks based on GGMs using an L1 graphical lasso for each chemical at each time point (i.e., 24 and 48 h) (Figures 2A, S4, and S5). The figures illustrate the estimated 190 partial correlation coefficients incorporated into the gene networks; edges with positive partial correlation values between two genes are shown in green, and edges with negative values are shown in red; the thickness, distance, and arrangement of the edges correspond to the degree of correlation between the two genes. Because these genes were obtained from the top five PCs that maximize the dispersion of the 24 chemicals using PCA, the estimated GGMs that describe the networks of all 20 genes differed considerably among chemicals, and it was difficult to classify the chemicals simply based on the network patterns as a whole. Therefore, for actual predictions, we decomposed the networks into their constituent edges rather than using them as a whole and used those with higher discriminative potential for training data as features. In other words, the partial correlation coefficients of the edges that are characteristic of the respective toxicity categories contributed to the SVM discrimination.

Prediction of six toxicity categories using KhES-3 cells and the GGM network
Using the 190 partial correlation coefficients in the GGM as input data, we predicted the toxicities of the 24 toxic chemicals in six categories using SVMs with leave-one-out-cross-validation (LOOCV). For the   iScience Article predictions, we followed the procedures described in a previous report on hEST-GN (Yamane et al., 2016) using four kernels (linear, polynomial, RBF, and maximum entropy) and increased the number of top features ranked by a t-test from 1 to 190. We also performed the prediction with the raw LFC values of the 3,200 (24 h) and 3,255 (48 h) transcription factor genes at each of the five concentrations relative to the vehicle-only expression. To compare the predictive accuracy, the same number of input data used for the GGM (i.e., up to 190 genes) was used as features. Predictions with the raw LFC values did not achieve a significantly higher prediction performance than the mean predictions using 10 uniform random numbers. On the other hand, in the prediction based on the GGM, the area under the receiver operating characteristic curve (AUC) values were R0.90 for chemicals in all toxicity categories, and because the prediction accuracy or AUC values was significantly high (p < 0.05), we concluded that prediction with high performance is possible. Predictions were performed separately at 24 and 48 h, but depending on the chemical, the time point at which higher prediction accuracy could be obtained differed; thus, neither time point was considered particularly superior in terms of yielding a better prediction. Overall, these results demonstrated that hEST-GN based on ES cell gene networks allows for the prediction of not only developmental toxicity but also broad toxicity categories including adult toxicity.
In addition, to compare with predictions based on the QSAR theory, we generated 5,666 molecular descriptors including 3D descriptors (Table S4) and performed predictions using top 1 to 190 feature genes according to the aforementioned method. None of the six categories, except for tubular nephrotoxin (accuracy of 91.7%), gave a significantly high prediction result. The results of all predictions are presented together in a table and as ROC curves (Table 2, Figures 2B and S6). These results suggest that chemical toxicity predictions that use the partial correlation coefficients of the GGM as features can achieve significantly higher accuracy than predictions based on gene expression values or QSAR. In the GGM-based prediction, the prediction accuracy for each of the 24 chemicals was examined from the SVM results. This analysis revealed that with respect to 16 chemicals (acetonylacetone, acrylamide, amitriptyline HCl, atorvastatin, bucillamine, chlorpheniramine, chlorpromazine, digoxin, doxorubicin, gentamicin, itraconazole, lithocholic acid, methapyrilene HCl, sunitinib, thioacetamide, and valproic acid), the prediction accuracy was 100% for all six categories at 24 and 48 h. On the other hand, for axitinib, cisplatin, and cyclosporin NT (13) HT (15) CT (13) GT (6) TT (7) NGC ( iScience Article A, the prediction accuracy was 66.6%, suggesting that the prediction of these chemicals is difficult (Table S5).
Pathway analysis of KhES-3 genes following exposure to 24 chemicals To determine the effects of exposure to chemicals on biological pathways, we performed a Hallmark pathway analysis by Gene Set Enrichment Analysis (GSEA) using all genes. When performing this analysis, we divided the samples into high dose (1/1 and ½ doses) and low dose ( 1 / 8 and 1/16 doses). For hepatotoxins, cardiotoxins, and glomerular nephrotoxins, differences were observed in the types of pathways that were induced or suppressed in comparison with other toxic chemicals in the high-dose samples at 24 h ( Figures 2C and S7). In addition, the responses of ES cell genes to the toxic chemicals were diverse and dependent on the type of chemical to which they were exposed and not limited to specific pathways such as apoptosis. This observation suggests that it is possible to predict toxicity categories on the basis of perturbed pathways that can be detected by transcription factors. On the other hand, differences in concentrations or among categories that may have been present at 48 h were not as pronounced as those at 24 h. However, analyses using available pathways based on human knowledge accumulated in the past provide limited information. Instead, the computational extraction of feature genes from the PCA of all genes without bias and predictions based on their GGM networks is likely to be more effective.

Selection of iPS cells as an alternative to human ES cells
The results of the present and previous studies suggest that hEST-GN can predict not only developmental toxicity but also broad toxicity categories with high performance. However, there are still hurdles to overcome, including ethical issues, before this system can be generally and widely accepted as a toxicity test. Accordingly, to make iPS cells a possible alternative to hEST-GN, we performed prescreening by comparing ATP assays with ES cells exposed to 20 toxic chemicals across a wide range of categories. As candidates, we used the top 20 cell lines selected from among Japanese male cell lines (Matsuda et al., 2020) derived from healthy individuals, which had been examined and ranked in terms of their differentiation potential into the three germ layers. For exposure concentration, we adopted the IC50 that was determined using the KhES-3 cell line and examined the toxicity response of human iPS cells. Among the candidate cell lines, we selected the top three cell lines with well-correlated growth rates at IC50 (HPS4138, HPS4234, and HPS4046) and confirmed the correlation coefficients of the growth rates at IC50 with KhES-3 using 20 of the 24 toxic chemicals investigated in this study. HPS4138 had the highest value of 0.94; accordingly, this cell line was used for the predictions as an alternative to ES cells (Table S6).

Prediction of chemicals in six toxicity categories using HPS4138 iPS cells
For HPS4138, which was selected by screening, we performed ATP assays with the 24 chemicals ( Figure S8; Tables S7 and S8). As in the case of ES cells, the cells were exposed to vehicle alone or five concentrations obtained by serial dilutions of stock solutions containing the maximum exposure concentration (i.e., the maximum value between IC0.1 and IC50 that did not cause an excessive degree of cell death); this experiment was repeated twice. Gene expression data for 20 genes selected from KhES-3 cells at 24 and 48 h were obtained by qRT-PCR, and GGMs were created for each of the 24 chemicals based on LFC values relative to the vehicle. Partial correlation coefficients were used for the prediction, as in the case of ES cells. In the prediction, data were created by integrating iPS cell data with ES cell data as well as by transductive transfer learning, in which toxicity category labels in the ES cell data were used for the learning to allow for category prediction using iPS cells. Assessment was performed by LOOCV, similarly to the predictions made using ES cells only. Chemicals in all categories except cardiotoxins and tubular nephrotoxins yielded AUC values from 0.82 to 0.99, and the accuracy or AUC was significantly higher than for results obtained with uniform random numbers. Thus, although this approach was not perfect, the results of the prediction using HPS4138 were very accurate for most toxicity categories (Table 3). The summary of toxicity category predictions for the 24 chemicals using HPS4138 is shown in Figure 3. Prediction was difficult for butylated HA, but satisfactory for the other chemicals (Table S9). These results demonstrate that if a gene expression database for toxicity responses could be created with ES cells, chemical toxicity prediction using iPS cells would also be possible by means of transductive transfer learning. Our findings also raise the possibility of achieving practical applications of toxicity testing using standardized or individualized iPS cells in the future. A description of transductive transfer learning is available on our web site at https://stempantox. stemcellinformatics.org.

DISCUSSION
Here, we presented a proof-of-concept study that enables a toxicity hazard assessment using human iPS cells and transfer learning based on the transcription factor gene network libraries made from the gene expression data of human ES cells exposed to 24 chemicals for 6 categories.
In this paper, we clarified that i) human ES cells are sufficient to detect not only developmental toxicities during embryogenesis but also broad toxicity categories such as adult toxins (neurotoxin, cardiotoxin, hepatotoxin, and glomerular and tubular nephrotoxins) and non-genotoxic carcinogens; ii) the chemical toxicity prediction using transcription factor gene networks of human ES cells shows an AUC = 0.90-1.00, which is significantly more accurate than predictions based on the QSAR theory or from raw gene expression data; iii) there exist differences in the biological pathways affected by the toxicity categories, suggesting the mechanisms that underlie the transcription factor networks that control the pathways may be used to predict toxicity categories; and iv) the gene network data from properly screened human iPS cells can successfully, although not perfectly, predict the toxicity categories at significant accuracies once the toxicities are learned by transfer learning using the models based on human ES cell data only.
Various alternative methods using pseudo-human systems, such as differentiated human cell lines (HepG2, MCF-7, HeLa, etc.), have been reported (Qu et al., 2021) since animal protection has become a higher priority in research. However, these lines are often derived from cancer or immortalized cells and thus have limited use (Kim et al., 2019). On the other hand, primary cells, which are assumed to resemble natural states in the human body, show batch-to-batch variability (Kim et al., 2019) and are difficult to collect at sufficient amounts. Furthermore, it is difficult to extrapolate toxicity tests on some cell types to other target cell types due to differences in cytotoxicity tolerance (Laschinski et al., 1991). Performing a multi-target toxicity prediction system based on stem cells, as we propose here, provides a more valid prediction of toxicity to a larger range of cell types.
Toxicological assessment using the transcriptome is frequently used by the U.S. EPA, Tox21 project and in Europe. Particularly, new approach methodologies (NAMs), which are any technology, methodology, approach, or combination thereof that can be used to provide information on chemical hazards and risk assessments that avoids the use of intact animals (New Approach Methods Work Plan: Reducing Use of Animals in Chemical Testing 2016), are often directional concepts using transcriptomics with other omics or traditional toxicology methods. Our study indicated that transcription factor gene networks exist in a master layer of biological pathways to activate molecular initiation events (MIEs), in which the initial chemical trigger starts an adverse outcome pathway (AOP) via DNA binding, receptor activation, or a disturbance of cellular/organelle systems (Allen et al., 2016), thus revealing toxicity reactions. Recent studies have iScience Article endorsed the idea that transcription factors of signal receptors might play a role in interfacing outside signals, such as aryl hydrocarbon or androgen, to activate toxicological AOPs in HepaRGä cells (Franzosa et al., 2021). In our system, we used 20 transcription factor genes from 5 PCs due to limited resources and cost, but it should be possible to customize the set of genes and PCs to reflect more accurately the specific AOPs in the endpoint organ. To pursue a full coverage of endpoint organs and AOPs, RNA-seq analysis and a library of all 4,033 transcription factor genes for the test chemicals at low cost are needed.
Previous systems using QSAR theory depend on the information of chemicals only and are thus inapplicable to mixtures such as food, Chinese or herbal medicines, and other compounds to assess toxicity as a whole. Our StemPanTox system detects the cellular toxic events of these mixtures, providing the iScience Article prediction of holistic cellular reactions, making it a resource for industries that mix independent chemical ingredients, including cosmetic, air-conditioner, and automotive companies who need to assess the toxicity of mixtures in their final or intermediate products, media, emissions, detergents, etc. In fact, more than 100 members from a wide variety of industry-government-academia fields are involved in our non-profit consortium (scChemRISC). By developing products from candidate substances that are predicted to have little toxicity, our system will contribute to industry not only for efficiency but also for human health.
Recently, toxicity reaction differences due to ethnicity, or genome haplotypes, have been widely reported due to the globalization of foods and products among countries. For example, catechin, which is contained in green tea and is widely consumed in Asian countries, is reported to induce severe liver injury in the United States (Oketch-Rabah et al., 2020). The CiRA Foundation (Kyoto, Japan) has announced myiPS cells, a project in which individuals can have their own iPS cells generated and banked.
Our StemPanTox system could allow a tailor-made chemical toxicity assessment for these cells to detect individual differences in toxic tolerance for different substances. Ideally, it also has the potential to reduce medical accidents if myiPS cells could be used to diagnosis whether a medication is toxic before receiving the treatment (Easley, 2019;Kim et al., 2019). In general, by performing a battery of toxicity assessments using multiple iPS cell lines from individuals with various haplotypes, our system may contribute to reducing toxicity accidents often caused by a small number of test samples of limited genomic variances.
In conclusion, the largest advantage of our StemPanTox is the ability to perform toxicity hazard assessments for multiple endpoints with high accuracy in a short amount of time and a low cost. We believe that our system will greatly benefit research aimed at amending gaps between current animal studies, which are resource-intensive and show uncertainty for human extrapolation, and NAMs, which are expected to use the latest science, technology, and computational models for improving regulatory toxicology guidelines.

Limitations of the study
Although we showed that human pluripotent stem cells are useful devices capable of qualitatively predicting adult toxicities of 24 chemicals in six categories, it is yet to be confirmed whether other chemicals and categories can be predicted and whether it is extendable to quantitative prediction.

STAR+METHODS
Detailed methods are provided in the online version of this paper and include the following:

DECLARATION OF INTERESTS
The authors declare no competing interests.  Since it has been reported that the toxicity of antioxidants such as catechin is suppressed in the presence of albumin (Zhang et al., 2016), maintenance culture was carried out for all cell lines including human ES cells using albumin-free Essential 8 Medium (Thermo Fisher Scientific) in six-well feeder-free culture dishes coated with 5 mg/mL vitronectin (VTN-N; Thermo Fisher Scientific). When seeding the cells, 10 mM CultureSureâ Y-27632 (FUJIFILM WAKO) was added, and medium exchange on day 1 and thereafter was performed without Y-27632.

Selection of toxic chemicals
Twenty-four chemicals were selected and mainly included neurotoxins, hepatotoxins, cardiotoxins, nephrotoxins, and non-genotoxic carcinogens ( Table 1). The presence or absence of toxicity was determined mainly on the basis of information regarding toxicity and assorted disorders and diseases available at PubChem (https://pubchem.ncbi.nlm.nih.gov/). With respect to neurotoxicity, among the chemicals previously reported in the literature, those having only developmental toxicity were classified as 'negative,' as the present study targeted adult toxicity. In addition, when determining the presence or absence of hepatotoxicity, chemicals with DILI rank R 3 were considered hepatotoxic chemicals; with regard to others, those with reliable reports of liver diseases were considered hepatotoxic. As for cardiotoxicity, chemicals that have been reported to be associated with heart disease were considered cardiotoxic. With regard to the kidney, due to its diverse and complex structure, the area of damage was divided into two sites, namely, the glomeruli and renal tubules. Chlorpheniramine and cyclopamine were the chemicals judged to be completely 'negative' and belonged to none of the toxicity categories examined in the present study.

Chemical exposures and determination of IC logistic model equations
DMSO or water was used as solvent (vehicle) for the 24 chemicals based on known information (https:// pubchem.ncbi.nlm.nih.gov/). For each of the 24 toxic chemicals, a stock solution was prepared with the highest soluble concentration. First, in order to perform the ATP assay to determine the exposure doses for testing, we performed 10 serial three-fold dilutions of the stock solution; the prepared exposure solution was added to the cells (i.e., exposure) so that the concentration in the medium was 0.1% of the exposure solution. The cells were cultured on 96-well black/clear flat bottom TC-treated plates (Falcon), 8,000 cells were seeded, and the medium was exchanged on day 1. The cells were exposed to chemicals on day 2. No medium exchange was performed after exposure, and the ATP assay was performed 48 h after exposure. For 100 mL of culture medium, 100 mL of CellTiter-Gloâ Luminescent Cell Viability Assay (Promega Corporation) was added, and emitted light was measured with a 2104 EnVision Multilabel Plate Reader (PerkinElmer). From four luminescence measurements for 10 concentrations and a blank, regression analysis was performed by fitting the three-parameter log-logistic model with the R-4.0.5 drc package, and IC0.1 and IC50 values were obtained (Tables S1 and S2). IC50 values were plotted using the ggplot2 package ( Figure 1B).

RNA-seq analysis
On day 2 after seeding, the cells were exposed to the chemicals. For each of the 24 chemicals, the exposure dose was set between IC0.1 and IC50 depending on the degree of cell death. iScience Article maximum exposure dose, five serial two-fold dilutions (1/1, 1/2, 1/4, 1/8, 1/16) were performed, and a total of six doses including a solvent-only control (vehicle) were used. After exposure, no medium exchange was performed, and samples were obtained at two time points (24 h and 48 h) with two repeats, i.e., a total of 24 x 6 x 2 x 2=576 samples. After RNA purification using an RNeasy Mini Kit (QIAGEN), sequence libraries were prepared for each sample using TruSeq Stranded mRNA Library Prep/TruSeq RNA Single Indexes Set A & Set B (Illumina, Inc.). For sequencing, high-throughput sequencing was performed using HiSeq4000 (Illumina, Inc.). We used bowtie-2.2.5 with the option "-very-sensitive-local" to map the obtained Illumina reads to Ensemble GRCh38r100 human cDNA and ncRNA sequences, added up the reads for each gene using MAPQ R 1 transcript, and obtained average counts of 28,652,809 and 28,291,682 reads for each sample at 24 h and 48 h, respectively. From these, we selected only transcription factor-related genes included in the Gene Ontology GO:0006351 using BioMart (4,032 genes). These genes were filtered using the statistical analysis language R-4.0.5 package edgeR (Robinson et al., 2010) (https://www.r-project.org/) with the filterByExpr function min.count = 30, min.total.count = 0, and then normalized to log2 counts per million (logCPM) using the voom function. Furthermore, the removeBatchEffect function was used to eliminate batch effects.

Differentially expressed gene analysis and principal component analysis
At 24 h and 48 h, for a total of 122 groups including 120 conditions (five concentrations each for 24 chemicals) and two solvent conditions (DMSO or water), we used a linear model fitting with the lmFit function of the limma package (Ritchie et al., 2015) in R and moderated t-statistics with eBayes to analyze differentially expressed genes (DEGs) with respect to gene expression levels in terms of the log-fold-change (LFC) between 24 chemicals and their corresponding solvent, and created a heatmap of genes with an LFC > 1 and FDR (false discovery rate) < 0.01 using the pheatmap package in R ( Figures 1C and S2). In addition, using the LFC values obtained for 120 conditions, PCA was performed for each of the two time points (24 h, 48 h) using the prcomp function in R ( Figures 1D and S3).

Gene network construction by Graphical Gaussian Model (GGM)
Based on results obtained in the aforementioned PCA, a total of 20 genes (two genes each with the top positive and negative loading values in the first to fifth PCs) were used to construct gene networks. To estimate the GGM for each of the 24 chemicals, we used the aforementioned LFC values to calculate the sparse partial correlation coefficient network with L1 graphical lasso using EBICglasso in the R package qgraph (https://cran.r-project.org/web/packages/qgraph/qgraph.pdf) (Epskamp and Fried, 2018). For model fitting, regular BIC with gamma = 0 was used, and regularization of sparsity was tried 1000 different ways with nlambda = 1000 for estimation. In the estimation, in order to avoid the problem of covariance matrices failing to be positive definite, calculations were performed with checkPD = FALSE (Figures 2A,  S4, and S5).

Prediction by support vector machine (SVM)
The SVM program and protocol used in the present study were adopted according to the report of Takahashi et al. (2018) Four kernel functions were used: linear, polynomial, RBF, and maximum entropy, and parameter types and combinations were calculated according to the above report (Takahashi et al., 2018). In the calculation, 190 values of partial correlation coefficients among 20 genes in the GGM for each of the 24 chemicals were used as input data, and using leave-one-out-cross-validation (LOOCV), the genes were ranked using the two-sample t-test (two-sided) in each iteration, and the maximum accuracy and the maximum AUC (area under the receiver operating characteristic curve) to achieve the maximum accuracy were recorded, varying the number of values from 1 to 190. For comparison, LFC data for the transcription factor genes at five concentrations before calculating the GGM (24 x 5 x 3,200 or 3,255 values in total) were used as input data, and the maximum accuracy and the corresponding AUC value were recorded in a similar manner. In addition, to compare with predictions based on QSAR, 5,666 molecular descriptors were created using alvaDesc (Affinity Science) (https://www.alvascience. com/alvadesc/). We obtained information from PubChem DB (https://pubchem.ncbi.nlm.nih.gov/) 3D Conformer data regarding 20 of the 24 chemicals; for the four other chemicals for which information could not be obtained from PubChem DB (cisplatin, cyclosporin A, digoxin, and gentamicin), the SMILES format was converted into 3D molecular descriptors using CORINA Classic (https://www.mn-am.com/ online_demos/corina_demo) and entered into alvaDesc.